rm(list = ls())
here::i_am(file.path("code", "06_iv_occupation.R"))
library(here)
source(here("code", "config.R"))

est_df <- read_parquet(here("data", "analysis", "analysis.parquet"))

occup_coefs <- map_dfr(
  unique(est_df$occ1950_c1930),
  function(x) {
    feols(
      outc ~ 1 |
        birthyr_cards + bpl_cards + board_identifier + exemption^married_cards +
        farmer + laborer + farm_laborer | vet_combined ~ order_num_serial_norm,
      data = est_df |>
        mutate(outc = occ1950_c1930 == x),
      cluster = ~ serial_num
    ) |>
      coeftable() |>
      as_tibble() |>
      mutate(occ = x)
  }
)

occ_xw <- read_csv(here("data", "xw", "occ1950.csv"))

p <- occup_coefs |>
  left_join(
    occ_xw,
    by = c("occ" = "value")
  ) |>
  filter(label != "N/A (blank)") |>
  filter(
    rank(Estimate) <= 25 | rank(-Estimate) <= 25
  ) |>
  mutate(
    bottom = Estimate < 0,
    label = gsub("\\(n.e.c.\\)", "", label)
  ) |>
  ggplot(aes(x = Estimate, xmin = Estimate - 1.96 * `Std. Error`,
             xmax = Estimate + 1.96 * `Std. Error`,
             y = fct_reorder(as.factor(label), Estimate))) +
  geom_vline(xintercept = 0, linetype = "dashed", alpha = .7) +
  
  my_theme() +
  theme(
    strip.text.x = element_blank(),
    legend.position = "none"
  ) +
  labs(x = "TSLS coefficient on veteran", y = "") +
  facet_wrap(~ bottom, ncol = 1, scales = "free_y")

p_color <- p + geom_point(aes(color = bottom), size = 2) + geom_linerange(aes(color = bottom))
ggsave(file.path(fig_dir, "color", "iv_occupation.pdf"), plot = p_color, width = 6, height = 9.5)

p_bw <- p + geom_point(size = 2, color = "black") + geom_linerange(color = "black")
ggsave(file.path(fig_dir, "bw", "iv_occupation.pdf"), plot = p_bw, width = 6, height = 9.5)
